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This  paper  presents  the  attitude  determination  method  for  the  Bifocal  Relay  Mirror 
Spacecraft  Simulator.  The  simulator  simulates  three-axis  motion  of  a  spacecraft  and  has  an 
optical  system  emulating  a  bifocal  space  telescope.  The  simulator  consists  of  three  control 
moment  gyroscopes,  rate  gyros,  two-axis  analog  sun  sensor,  and  two  inclinometers.  The 
five-foot  diameter  platform  is  supported  on  a  spherical  air  bearing  to  offer  a  low-torque 
environment.  This  paper  demonstrates  two  attitude  determination  methods  employing 
the  measurements  from  a  two-axis  analog  IR  sensor,  two  inclinometers,  and  a  triaxial 
gyroscope.  The  first  method  implements  the  conventional  Kalman  filter  algorithm.  The 
second  method  uses  a  nonlinear  observer  derived  from  the  Lyapunov’s  direct  method. 
Analytical  and  experimental  results  are  presented  to  validate  the  proposed  algorithm. 

I.  Introduction 

The  Naval  Postgraduate  School  (NPS)  Three- Axis  Spacecraft  Simulator  (TASS),  which  is  shown  in  Fig.  1, 
is  designed  for  developing  and  validating  acquisition,  tracking,  and  pointing  technologies  for  Bifocal  Relay 
Mirror  Spacecraft  (BRMS).  The  BRMS  consists  of  single  axis  gimballed  receive  and  transmit  telescopes.  The 
telescopes  are  used  to  redirect  the  laser  light  from  the  ground  (or  space)  to  a  distant  target  which  necessitate 
a  tight  pointing  requirement.1 

The  TASS  simulates  spacecraft  three-axis  motion  and  has  an  optical  system  simulating  a  single  space 
telescope.  The  five-  foot  diameter  TASS  platform  is  supported  on  a  spherical  air  bearing  to  offer  a  low- 
torque  environment.  The  simulator  consists  of  three  variable  speed  control  moment  gyroscopes  as  actuators, 
rate  gyros  for  angular  rate,  two-axis  analog  sun  sensor  and  two  inclinometers  for  attitude  sensors.  Detailed 
description  of  the  TASS  hardware  is  given  in  Ref.  2. 

This  paper  focuses  on  the  development  and  validation  of  the  attitude  determination  algorithm  of  the 
TASS  using  two  inclinometers,  two-axis  IR  sensor,  and  a  triaxial  rate  gyroscope.  Although  inclinometers 
cannot  be  used  in  a  spacecraft,  they  can  still  be  employed  in  testbeds  to  get  useful  attitude  information  with 
good  accuracy  at  low  cost.  Also,  if  the  light  source  is  not  far  enough,  the  effect  of  transitional  motion  of  the 
IR  sensor  during  the  rotation  of  the  testbed  should  be  taken  into  account  when  calibrating  the  IR  sensor. 

Quaternion  attitude  parametrization  offers  a  singularity-free  method  to  represent  three-dimensional  rota¬ 
tions  using  four-dimensional  unit  vectors.  But  the  four  elements  of  the  unit  quaternion  cannot  be  estimated 
independently  because  of  the  normality  constraint. 

It  is  a  standard  approach  to  include  the  gyro  bias  in  the  overall  kinematics  model  and  to  estimate  the 
actual  attitude  by  Extended  Kalman  Filter  (EKF)  techniques.  The  EKF,  as  an  extension  of  the  Kalman 
Filter  (KF)  to  nonlinear  models,  works  satisfactorily  when  the  linearized  model  exhibits  sufficient  accuracy. 
Since  the  EKF  utilizes  a  linearization  at  each  step,  the  nonlinear  normality  constraint  of  the  quaternions 
cannot  be  rigidly  enforced  within  the  EKF  computational  structure.  Various  methods  have  been  developed 
to  overcome  this  deficiency  and  are  described  in  Refs.  3,  4,  5,  and  6. 

The  nominal  attitude  determination  method  used  for  the  TASS  is  a  basic  six-state  attitude  estimation 
filter  known  as  the  Multiplicative  EKF  (MEKF).4  In  the  MEKF  formulation  of  the  TASS,  the  Modified  Ro¬ 
drigues  Parameters  (MRPs)  are  used  as  an  unconstrained  representation  of  the  attitude  error  with  quaternion 
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Figure  1.  NPS  Three-Axis  Spacecraft  Simulator  (TASS) 


representing  the  reference  attitude  of  the  defined  attitude  error.  The  main  advantage  of  using  MRPs  over 
other  three  dimensional  attitude  parameterizations  is  that  it  has  a  singularity  at  360°. 

In  this  research,  we  propose  a  different  approach  to  estimate  the  attitude  and  gyro  bias.  The  proposed 
method,  based  on  a  Lyapunov  observer,  employs  both  the  quaternion  and  the  MRPs,  and  it  fully  exploits  the 
nonlinear  structure  of  the  attitude  kinematics.  Consequently,  global  convergence  is  analytically  guaranteed, 
which  means  that  the  system  can  recover  from  large  errors. 

In  what  follows,  the  sensor  models  are  given,  then  the  MEKF  is  reviewed.  Next,  a  nonlinear  Lyapunov 
observer  is  derived.  Finally  numerical  and  experimental  validation  of  the  proposed  algorithm  is  presented. 


II.  Sensors  Model 

For  the  TASS  attitude  determination,  two  inclinometers  are  mounted  orthogonally  to  each  other  in  the 
testbed’s  body  x  and  z  axis  as  shown  in  Fig.  2  to  measure  the  direction  of  the  gravity  vector.  The  two-axis 
IR  sensor  is  placed  near  the  edge  of  the  testbed  parallel  to  the  body  x  axis  vector.  Since  the  IR  sensor  is  not 
situated  at  the  center  of  the  rotation  point,  the  effect  of  the  transitional  motion  of  the  sensor  with  respect 
to  the  light  source  may  have  significant  effect  on  the  sensor  reading. 

In  this  section,  we  discuss  in  detail  how  we  can  extract  gravity  vector  with  two  inclinometers  then  explain 
a  simple  way  to  compensate  the  transitional  motion  of  the  IR  sensor  on  the  sensor  reading. 


A.  Inclinometers 


The  inclinometers  mounted  on  the  TASS  use  capacitive  liquid  to  measure  the  tilt  angle  of  the  platform 
with  respect  to  the  gravity  vector.  The  measuring  range  of  the  sensors  is  ±  30°  with  a  linear  angle  output 
accuracy  of  0.2°.  The  transverse  axis  sensitivity  is  about  1%  of  the  angle  reading  at  30°  tilt.7  Figure  3  shows 
a  diagram  of  an  inclinometer  fixed  in  the  testbed  that  is  tilted  by  +6bi-  The  vector  Li  lies  along  the  liquid 
surface  and  it  is  perpendicular  to  the  gravity  vector.  The  right-handed  sensor  coordinate  frame  i\  is  defined 
by  x/i,  y/i,  and  z/i.  The  direction  of  the  vector  Li  can  be  expressed  in  the  sensor  frame  I\  as 


JlL 


1 


—  tan(0n) 


(1) 
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Figure  2.  Testbed  placement  in  the  NPS  lab 


with  the  assumption  that  its  Xu  component  equals  to  one.  Using  a  similar  approach  for  the  second  incli¬ 
nometer,  the  unit  gravity  vector  can  be  found  by  taking  the  cross  product  between  Li  and  L2  as 

Bg=BL!  xbL2/|L!  xL2|  (2) 

after  converting  and  L2  from  the  sensor  frames  (. I\  and  /2)  into  the  testbed’s  body  frame  B.  With  small 
angle  approximations,  the  Li  and  L2  output  angles  correspond  to  roll  and  pitch,  respectively. 

B.  IR  Sensor 

The  IR  sensor  is  mounted  near  the  edge  of  the  testbed  and  senses  the  20  W  lamp  located  at  110  cm  above 
the  floor  and  3  m  away  from  the  center  (Ob)  of  the  testbed  body  frame.  A  top  view  of  the  IR  sensor  and 
the  incident  light  with  the  testbed  is  shown  in  Figs.  4.  From  figure  4,  we  can  see  that  the  translational 
displacement  of  the  IR  sensor  cannot  be  neglected  since  s  and  s'  are  not  parallel  with  each  other.  The 
relation  between  these  two  vectors  can  be  expressed  as 

s'  =  s  +  rIR  (3) 

where  s  is  the  vector  from  the  sensor  to  the  source  light  and  r jR  is  the  position  vector  of  the  IR  sensor. 
The  vector  from  the  testbed  center  of  rotation  to  the  source  light  is  expressed  as  s'.  For  the  attitude 
determination,  the  physical  quantity  we  need  as  an  observation  vector  is  the  unity  vector  of  s'  instead  of  the 
measured  vector  s.  Equation  3  can  be  rewritten  as 

s'  =  Ls  +  rIR  (4) 

where  L  is  a  scalar  representing  the  magnitude  of  the  vector  s  and  can  be  calculated  by  taking  the  norm  on 
both  sides  of  Eq.  4  giving 

|s'|  =  |Ls  +  r/fl|  (5) 
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Inclinometer 


Figure  3.  Inclinometer  single-axis  view 


where  the  quantities  s,  iq^  and  |s'|  are  known.  Then  letting  the  known  measured  vector  s  =  [x,  y ,  z]T  and 
the  position  of  the  IR  sensor  r m  =  [a,  &,  c]T,  the  unknown  magnitude  L  can  be  solved  as 

zc  +  yb  +  xa  —  \j2zcyb  +  2  zcxa  +  2  ybxa  +  z2(—b2  —  a2  +  d2)  +  y2{—c2  —  a2  +  d2)  +  x2(—c2  —  b2  +  d 2) 


z2  +  y2  +  x2 

Therefore  the  unit  vector  s'  can  be  calculated  using  Eqs.  5  and  6. 


(6) 


III.  Attitude  Estimation  Using  Extended  Kalman  Filter 

A  widely  used  six-state  attitude  estimation  method  known  as  the  MEKF4  is  implemented  for  the  attitude 
determination  of  the  TASS.  In  the  MEKF  of  the  TASS,  Modified  Rodrigues  Parameters  (MRPs)  are  used 
as  an  unconstrained  representation  of  the  attitude  error,  with  some  quaternion  representing  the  reference 
attitude  of  the  defined  attitude  error.  The  state  vector  of  the  MEKF  consists  of  six  elements,  three  for  the 
MRPs  and  three  for  the  gyro  bias. 

A  quaternion  is  defined  as 

qi3  {t) 

94  (t) 


q  (t)  = 


(7) 


where  the  vector  q13  is 


and  the  scalar  q±(t)  is 


qi3(i)  = 


Qi(t) 

92  {t) 

93  W 


=  n  sm 


(f 


.  ( 9(t) 

q4{t)  =  COS  I  — 


(8) 


(9) 


where  n  is  a  unit  vector  indicating  the  principal  rotation  axis  and  9  is  the  principal  rotation  angle.  The 
quaternion  components  satisfy  the  following  normalization  constraint: 


T  2  ,  2  ,  2  ,  2  i 

q  q  =  9i  +  92  +  93  +  94  =  1 


A  three  dimensional  MRPs  vector  is  defined  as 


»«)=  4q'a(i) 


1  +  94(0 

In  the  MEKF  formulation,  the  true  attitude  quaternion  q  is  represented  by 

q (t)  =  tfq(a(t))  <8>  qref(t) 


(10) 


(11) 


(12) 
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where  q ref  and  <5q(a(f))  are  the  reference  and  the  error  quaternion  in  the  body  frame,  respectively.  The 
quaternion  multiplication  ®  is  defined  as4 


q'(g>q  = 


<3^13  +  94qi3  -  q'i3  x  qi3 

9494  -  qi3  ■  qi3 


(13) 


The  error  quaternion  (5q(a(t)),  which  represents  the  error  between  the  true  and  the  reference  attitude,  can 
be  expressed  in  terms  of  MRPs  as4 


Sq{a{t)) 


1 

16  +  a(t )2 


8a(t)  \ 

16  —  a[t)2  J 


where  a  is  the  magnitude  of  the  MRPs  vector. 

The  quaternion  kinematic  equations  of  motion  are  given  by 


q(*)  = 


1 

2 


uj(t) 

0 


®  q (t) 


(14) 


(15) 


where  u>  is  the  angular  velocity  of  the  body  frame.  In  the  MEKF,  a  is  chosen  to  keep  a  =  0  during  the 
state  propagation  so  that  the  quaternion  q ref  represents  an  estimate  of  correctly  normalized  true  attitude 
estimate. 

The  gyroscope  model  is  given  by 


uj(t)  =  u>(t)  -  b{t)  -  r]v(t)  (16) 

b(t)  =  Vn{i)  (17) 

where  ih(t)  represents  the  gyro  output,  b(t)  is  the  gyro  bias,  f)v(t)  and  are  uncorrelated  zero-mean 

white  noise  processes. 

The  Kalman  filter  uses  quaternion  kinematics  Eq.  15  with  q  =  qre/  and 

u i{t)  =  d>(t)  =  u)(t)  —  b(£)  (18) 

b  (t)  =  0  (19) 
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to  propagate  the  filter  state. 

The  error  covariance  propagation  matrix  V  satisfies 

V{t)  =  F  [x(t),  t]  V(t)  +  V(t)F  [x(t),  tf  +  Q 


where  F(t)  and  G(t)  matrices  are  given  by 

m  = 

G(t)  = 


-[w(t)x]  -J3x  3 
03x3  03x3 

^3x3  03x3 

03x3  hx3 


and  Q(t)  is  the  process  noise  covariance  matrix  of  the  following  form 


Q 


(20) 

(21) 

(22) 

(23) 

(24) 


where  the  terms  of,  and  of  are  the  variances  of  //„  and  77^,  respectively. 

The  measurement  update  is  done  with  the  two  vector  measurements  from  the  IR  sensor  and  the  incli¬ 
nometers.  The  sensitivity  matrix  is  given  by 


H  = 


[^(q)Jvgx] 

[A(q)^s'x] 


03x3 

03x3 


q=qre/ 


(25) 


where  A  is  the  attitude  matrix  of  the  body  frame  with  respect  to  an  inertial  frame.  The  vectors  wg  and  Ns' 
are  the  inertial  frame  representations  of  the  vectors  g  and  s'. 

Bad  measurement  data  (large  spikes),  caused  by  disturbances  from  the  testbed  electrical  components, 
are  removed  by  monitoring  the  filter  innovation  vector  residual  before  they  are  processed  by  the  filter. 

For  further  details  on  the  MEKF,  refer  to  Refs.  3,  4,  and  6. 


IV.  Nonlinear  Gyroscope  Bias  Observer 

In  this  section,  a  nonlinear  attitude  and  bias  observer  is  derived  based  on  Lyapunov’s  direct  method.  Like 
the  MEKF,  this  observer  estimates  attitude  and  gyro  bias  by  employing  both  the  quaternion  and  the  MRPs. 
It  fully  exploits  the  nonlinear  structure  of  the  attitude  kinematics.  As  a  consequence  global  convergence  is 
guaranteed  at  least  analytically,  which  means  that  the  system  can  recover  from  large  errors. 

To  design  the  nonlinear  observer  consider  Eq.  12.  The  reference  quaternion  q ref  is  defined  such  that  it 
follows 

®q  refit)  (26) 

where  the  reference  body  rate  is  now  defined  as 

uref{t)  =  G){t)-b(t)-e{t)  (27) 

where  e  is  a  3  by  1  vector,  which  will  be  defined  later.  The  MRPs  (a)  satisfy  the  following  kinematic 
equations:4 

a  =  [(1  -  —  a2)I3x3  +  -aaT](w  -  u >ref)  +  -a  x  (w  +  uref)  (28) 

where  u>  is  the  true  body  rate. 

Consider  the  Lyapunov  function  candidate  of  the  form 

V(a,  b)  =  ^aaTa+  ^/3bT b  (29) 


q  refit)  =  2 


Wref{t) 

0 
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where  a  and  [3  are  positive  constants,  and  the  gyro  bias  estimate  error  b  =  b  —  b.  Then  the  time  derivative 
of  the  Lyapunov  function  V  is  given  by 

y(a,  b)  =  a-aTa  +  (dbTb 

=  aaT[(l  -  ~^a2)I3x3  +  ^aaT](u;  -  ujref)  +  /3bTb  (30) 

where  Eq.  28  is  substituted  and  aT[a  x  (u>  +  o>re/)]  =  0  is  used.  Using  Eqs.  16  and  27,  and  neglecting  the 
noise  rjv,  the  vector  (a;  —  u)ref)  in  Eq.  30  can  be  expressed  as  u)  —  c oref  =  b  +  e.  Then  we  can  rewrite  Eq.  30 
as 

U(a,  b)  =  a-aT[(l  —  Y^a2)I3x3  +  ^aaT](b  +  e)  +  /3bTb 

=  a[l  +  ^a2]aTb  +  a[l  +  ^a2]aTe  + /3bTb  (31) 

If  we  choose 

e  =  -Aa  (32) 

~  (V  1 

b  =  -^[1+^a  (33) 

where  A  is  a  positive  constant,  then  V  can  be  simply  expressed  as 

V  =  — Aa[l  +  ^a2]a2  (34) 

which  is  negative  semi-definite. 

Because  V  <  0,  the  observer  is  only  stable  in  the  sense  of  Lyapunov  but  not  asymptotically  stable. 
To  prove  the  asymptotic  stability  of  the  observer  about  a  =  0  and  b  =  0,  we  use  theorem  8.5  of  Ref.  8. 
According  to  the  theorem,  asymptotic  stability  can  be  achieved  1)  if  the  first  k  —  1  derivatives  of  U(x)  is 
zero  when  evaluated  on  the  set  Cl,  which  is  the  set  of  non-empty  state  vector  satisfying  x  £  =>  U(x)  =  0, 

and  2)  the  kth  derivative  of  V (x)  is  negative  definite  on  the  set  Cl  with  an  odd  number  k. 

Therefore  the  set  Cl  can  be  calculated  from  Eq.  34  and  is  given  by  Cl  =  {(a,  b)|a  =  0}.  The  second 
derivative  of  V  is 

V  =  —  Aa[2  +  ^a2]aT[(l  -  ^ a2)I3x3  +  ^aaT](b  +  e)  (35) 


which  is  still  0  when  evaluated  at  fl  =  {(a,  b)|a  =  0}.  Taking  the  time  derivative  of  Eq.  35,  the  third 
derivative  of  V  can  be  written  as 


V  =  -AaiaTaaT[(l  -  ^ a2)I3x3  +  ^aaT](b  +  e)  -  Aa(2  +  ia2)aT[(l  -  ^ a2)I3x3  +  ^aaT](b  +  e) 

-  Aa(2  +  |fl2)aT(-^aTa /  +  ^aaT  +  ^aaT)(b  +  e)  -  Aa[2  +  ^a2] aT[(l  -  ^a2)I3x3  +  ^aaT](b  +  e) 

(36) 


Then  evaluating  Eq.  36  on  the  set  Cl  and  simplifying  gives 

— 2Aa(b  +  e)T(b  +  e)  (37) 

which  is  negative  definite  for  all  fi  =  {(a,  b)|a  =  0}.  Therefore  we  can  conclude  that  Eqs.  28  and  33  are 

globally  asymptotically  stable  with  no  noise.  Assuming  that  the  bias  is  nearly  constant  then  b  ~  b  can  be 
used  to  estimate  the  bias  with  Eq.  33.  Note  that  the  attitude  errors  represented  by  MRPs  can  be  obtained 
from  Eq.  12  with  the  assumption  that  q  ~  q  where  q  is  the  quaternion  measurement  converted  from  the 
two  vector  observations.  In  this  paper,  we  use  the  deterministic  TRIAD  algorithm9  to  determine  an  attitude 
from  the  gravity  and  the  sun  vector.  In  the  TRIAD,  the  unit  gravity  vector  is  used  as  the  first  basis  since 
it  is  more  accurate  then  the  sun  vector  which  is  from  the  IR  sensor.  Note  that  when  noises  are  present  in 
the  measurements,  estimated  values  still  converge  to  the  region  around  the  true  values  within  error-bounds 
which  are  related  to  the  noise  upper  bounds. 
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V.  Simulation  and  Experiment  Results 


In  this  section,  several  numerical  simulations  and  experiments  have  been  done  to  validate  the  MEKF  and 
the  Lyapunov  observer  in  the  presence  of  noise  in  the  sensors  data.  The  TASS  has  a  Northrop  Grumman 
Litton  LN-200  IMU  consisting  of  three  fiber  optics  rate  gyroscopes  and  a  coarse  two-axis  analog  sun  sensor 
consisting  of  four  photo-transistors. 

A  20  W  halogen  source  light  is  placed  on  the  wall  3  m  away  from  the  TASS  center  of  rotation.  The 
gyro  measurement  data  are  simulated  using  Eqs.  16  and  17  with  the  noise  standard  deviations  of  <Jupsuon  = 
0.0195 ° /s1/2  and  =  0.0057 °/s3/2.  Initial  gyro  bias  of  b(t0)  =  [4.3  3.2  2.3]T  *  10”4  (rad/s)  is  used.  For 
the  simulations,  a  Gaussian  white  noise  of  cr  =  0.001  has  been  added  to  the  measurement  of  the  gravity  and 
the  sun  vector.  The  initial  estimate  for  the  attitude  quaternion  has  been  set  to  q  =  [0  0  1  0]T  (whereas  the 
true  attitude  is  q  =  [0  0  0  1]T)  and  the  bias  has  been  set  to  b  =  [0  0  0]T.  Data  sampling  interval  of  0.025  s 
has  been  chosen  for  all  sensors  with  a  4th  order  Runge-Kutta  method  for  the  numerical  integration  method. 
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Figure  7.  Effect  of  ^  changes  on  the  Lyapunov  observer  estimate 


Figures  5  and  6  show  the  comparisons  of  bias  and  attitude  estimate  computed  using  the  MEKF  and  the 
nonlinear  Lyapunov  observer.  For  the  Lyapunov  observer,  A  and  of  Eqs.  32  and  33  are  chosen  to  be  0.1 
and  0.005,  respectively.  We  can  see  from  the  figures  that  the  Lyapunov  observer  estimate  accurately  follows 
the  estimate  of  the  MEKF  after  200  s. 

The  effects  of  changes  in  the  values  of  ^  and  A  on  the  estimation  performance  have  been  shown  in  Figs.  7 
and  8,  respectively.  We  can  see  from  Fig.  7  that  as  ^  increases  the  bias  estimate  becomes  more  sensitive  to 
the  sensor  noise.  From  Fig.  8,  we  can  see  that  the  value  A  should  be  kept  small  to  get  a  good  performance 
in  the  presence  of  sensor  noise. 

For  the  experiments,  standard  deviations  of  the  Gaussian  white  noises  cr  =  [0.002  2  x  10~6  0.001] r 
(gravity  vector)  and  cr  =  [0.001  0.0099  0.0358]T  (sun  vector)  are  used  to  tune  the  MEKF.  Figure  9  shows  the 
result  of  the  bias  estimate  using  both  the  MEKF  and  the  Lyapunov  observer.  We  can  see  from  this  figure 
that  the  estimated  values  for  both  algorithms  are  consistent. 
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VI.  Conclusion 


In  this  paper  the  attitude  determination  method  for  the  NPS  Bifocal  Relay  Mirror  Spacecraft  Simula¬ 
tor  has  been  shown.  The  simulator  simulates  three-axis  motion  of  the  spacecraft.  The  nominal  attitude 
determination  approach  is  the  MEKF  based  on  the  modified  Rodrigues  parameters  (MRPs)  as  an  uncon¬ 
strained  representation  of  the  attitude  error.  The  MEKF  uses  measurements  from  the  two-axis  analog  IR 
sensor,  two  inclinometers,  and  gyroscopes.  A  method  to  extract  the  gravity  vector  from  inclinometers  has 
been  shown  with  a  compensation  method  for  the  effect  of  transitional  motion  of  the  IR  sensor  due  to  the 
rotation  of  the  testbed.  As  an  alternative  for  the  MEKF  method,  a  Lyapunov  based  nonlinear  observer 
has  been  derived  which  provides  a  globally  asymptotically  stable  algorithm.  The  nonlinear  observer  works 
well  even  in  the  presence  of  noises  and  large  initial  errors.  In  addition,  it  is  easy  to  operate  since  we  have 
only  a  few  observer  parameters  to  tune  compared  to  the  MEKF.  Moreover  the  algorithm  is  very  simple  and 
computationally  inexpensive.  Numerical  and  experimental  validations  have  been  shown  for  both  attitude 
determination  algorithms.  The  main  drawback  of  the  nonlinear  observer  is  that  it  needs  high  data  sampling 
frequency. 
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